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Une methode semi-Lagrangienne en avant pour la resolution 
numerique de 1 'equation de Vlasov 

Resume : Ce document concerne la resolution numerique de I'equation de Vlasov qui est un modele cinetique 
permettant de decrire devolution d'un plasma. Elle est couplee a une equation, par exemple I'equation de 
Poisson, permettant le calcul des champs auto-consistants. Le modele couple est non lineaire. Nous proposons 
ici une nouvelle methode semi-Lagrangienne dans laquelle les caracteristiques sont integrees en avant, dans 
le sens du temps, contrairement a la methode semi-Lagrangienne classique qui integre les caracteristiques en 
arriere. L'autre ingredient de la methode semi-Lagrangienne est une technique dc deposition de I'information 
portee par les particules sur une grille. Celle-ci est basee sur un produit tensoriel de splines cubiques de 
sorte qu'en deux dimensions la formule de deposition implique 16 points de la grille. Grace a cette nouvelle 
technique la methode semi-Lagrangienne devient completement explicite ct pcut s'appuyer sur des solveurs 
numeriques d'equations differentielles classiques, rendant ainsi plus simple la montee en ordre en temps. La 
methode est validee sur plusieurs cas tests representatifs des problemes realistes bases sur ce modele. 

Mots-cles : methode semi-Lagrangienne, Runge-Kutta, simulation de plasmas, equation de Vlasov 
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1 Introduction 

Understanding the dynamics of charged particles in a plasma is of great importance for a large variety of 
physical phenomena, such as the confinement of strongly magnetized plasmas, or laser-plasma interaction 
problems for example. Thanks to recent developments in computational science and in numerical methods, 
meaningful comparisons between experience and numerics are becoming possible. 

An accurate model for the motion of charged particles, is given by the Vlasov equation. It is based on 
a phase space description so that non-equilibrium dynamics can accurately be investigated. The unknown 
f{t,x,v) depends on the time t, the space x and the velocity v. The electromagnetic fields are computed 
self-consistently through the Maxwell or Poisson equations, which leads to the nonlinear Vlasov-Maxwell or 
Vlasov-Poisson system. 

The numerical solution of such systems is most of the time performed using Particle In Cell (PIC) 
methods, in which the plasma is approximated by macro-particles (see [3]). They are advanced in time 
with the electromagnetic fields which are computed on a grid. However, despite their capability to treat 
complex problems, PIC methods are inherently noisy, which becomes problematic when low density or highly 
turbulent regions are studied. Hence, numerical methods which discretize the Vlasov equation on a grid of 
the phase space can offer a good alternative to PIC methods (see [SJ [5J IHl UHl H] ) ■ The so-called Eulerian 
methods can deal with strongly nonlinear processes without additional complexity, and are well suited for 
parallel computation (see [12] )■ Moreover, semi-Lagrangian methods which have first been introduced in 
meteorology (see [13 [211 HI] ) , try to take advantage of both Lagrangian and Eulerian approaches. Indeed, 
they allow a relatively accurate description of the phase space using a fixed mesh and avoid traditional 
step size restriction using the invariance of the distribution function along the trajectories. Standard semi- 
Lagrangian methods calculate departure points of the characteristics ending at the grid point backward in 
time; an interpolation step enables to update the unknown. 

In this work, we consider the numerical resolution of the two-dimensional Vlasov equation on a mesh of the 
phase space using a forward semi-Lagrangian numerical scheme. In the present method, the characteristics 
curves are advanced in time and a deposition procedure on the phase space grid, similar to the procedure 
used in PIC methods for the configuration space only, enables to update the distribution function. 

One of the main cause for concern with semi-Lagrangian methods, is computational costs. With Back- 
ward semi-Lagrangian methods (BSL), the fields have to be computed iteratively, with Newton fixed point 
methods, or prediction correction algorithms. This is due to an implicit way of solving the characteristics 
(see fT9\ for details). This strategy makes high order resolution quite difficult and expensive. Making the 
problem explicit enables to get rid of iterative methods for the characteristics, and to use for example high 
order Runge Kutta methods more easily. This is one of the main advantages of the present Forward semi- 
Lagrangian (FSL) method. Once the new position of the particles computed, a remapping (or a deposition) 
step has to be performed. This issue is achieved using cubic spline polynomials which deposit the contri- 
bution of the Lagrangian particles on the uniform Eulerian mesh. This step is similar to the deposition 
step which occurs in PIC codes but in our case, the deposition is performed in all the phase space grid. 
Similarities can also be found in strategies developed in [HI [TH [THj for meteorology applications. 

In order to take benefit from the advantages of PIC and semi-Lagrangian methods, and since the two 
methods (PIC and FSL) really look like each other, except the deposition step, we have also developed 
a hybrid method, where the deposition step is not performed at each time step, but every T time steps. 
During the other time steps, the fields are computed directly at the new position of the particles. It shall 
be noticed that however the present method is not a real PIC method, since the particle weights are not 
constant. Indeed, in this method based on a description of the unknown using cubic spline polynomials, 
the spline coefficients play the role of the particle weights, and are updated at each time step. This kind of 
hybrid approach has been developed recently in a slightly different framework in ^20] inspired by [7] . 

This paper is organized as follows. In the next section, the two Vlasov equations which will be dealt with 
are presented. Then, we shall introduce the Forward semi-Lagrangian (FSL) method, always regarding it 
comparatively to Backward semi-Lagrangian (BSL) methods. Afterward, numerical results for several test 
cases are shown and discussed. Eventually, some specific details are given in two appendices, one for the 
computation of an exact solution to the Landau damping problem, and the other about the solution of the 
Poisson equation for the Guiding-Center model. 
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2 Models in plasma physics 

In this section, we briefly present two typical reduced models from plasma physics for the description of the 
time evolution of charged particles. These two-dimensional models are relevant for more complex problems 
we are interested in and shall be used to validate our new method. 



2.1 Vlasov-Poisson model 



We consider here the classical ID Vlasov-Poisson model, the unknown of which / = f{t, x, v) is the electron 
distribution function. It depends on the space variable x £ [0, L] where L > is the size of the domain, 
the velocity variable u G E. and the time t > 0. The Vlasov equation which translates the invariance of the 
distribution function along the characteristics then writes 



Of 

-^+vd,f + E{t,x)d,f^O, 



(2.1) 



with a given initial condition f{0,x,v) 
thanks to the distribution function / 



fo{x,v). The self-consistent electric field E{t,x) is computed 



d,E{t, x) 



f{t,x,v)dv - Pi 



R 



E{t,x)dx = 0, 



(2.2) 



where pi denotes the ion density which forms a uniform and motionless background in the plasma. 

The Vlasov-Poisson model constitutes a nonlinear self-consistent system as the electric field determines 
/ with (|2.ip and is in turn determined by it in (|2.2p . It presents several conserved quantities as the total 
number of particles, the norms {p > 1) defined by ||/||lp = iff l/l^dxdw)^/^, the momentum and the total 
energy, as follows: 



dec d 

- // f{t,x,v)dxdv ^ -\\f{t)\\LP 



d_ 
di 
d_ 
di 
0. 



vf(t, X, v)dxdv 
v^fit,x,v)dxdv + lEit,xrdxdv 



One of the main features of the present work is to develop accurate numerical methods which are able to 
preserve these exactly or approximately these conserved quantities for long times. 



2.2 Guiding-center model 

We are also interested in other kinds of Vlasov equations. For instance, in the guiding-center approximation. 
Charged particles in magnetized tokamak plasmas can be modeled by the density / = f{t,x,y) in the 2 
dimensional poloidal plane by 

^ + E^{x,y)-\/f = 0, (2.3) 

coupled self-consistently to Poisson's equation for the electric field which derives from a potential $ — ^{x,y) 

- A$(i, X, y) = fit, X, y), E{t, x, y) = -V$(t, x, y). (2.4) 

In equation (12. 3p . the advection term E^ — {Ey,—Ex) depends on {x,y) and the time-splitting cannot be 
simply applied like in the Vlasov-Poisson case. Hence, this simple model appears to be interesting in order 
to test numerical methods. 

The guiding-center model (|2.3p - l|2.4p also presents conserved quantities as the total number of particles 
and norm of / (energy) and E (enstrophy) 
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2.3 Characteristic curves 

We can re-write Vlasov equations in a more general context by introducing the characteristic curves 

^ = UiXit),t). (2.6) 

Let us introduce X(t,x,s) as the solution of this dynamical system, at time t whose value is x at time s. 
These are called the characteristics of the equation. With X{t) a solution of (|2.6p . we obtain: 

J,ifim,t)) - I + f • V./ = I + UiXitU) . V./ = 0. (2.7) 

which means that / is constant along the characteristics. Using these notations, it can be written 

fiX{t; X, s), t) = f{X{s; x, s), s) = f{x, s) 

for any times t and s, and any phase space coordinate x. This is the key property used to define semi- 
Lagrangian methods for the solution of a discrete problem. 



3 The forward semi-Lagrangian method 

In this section, we present the different stages of the forward semi-Lagrangian method (FSL) and try to 
emphasize the differences with the traditional backward semi-Lagrangian method (BSL). 



3.1 General algorithm 

Let us consider a grid of the studied space (possibly phase-space) with N^: and Ny the number of points in 
the X direction [0, L^] and in the y direction [0, Ly]. We then define 

Ax = La;/N.^, Ay = Ly/Ny, Xi = iAx, y^ = jAy, 

for i — 0,..,Nx and j — 0,..,Ny. One important point of the present method is the definition of the 
approximate distribution functions which are projected on a cubic B-splines basis: 

fit, x,y)^Y. ^''^(^ - Xk,yun)S{v - X^it; Xk, yi,n), (3.8) 

k,l 

where X{t; Xk,yi,t^) — {Xi,X2){t,Xk,yi,t") corresponds to the solution of the characteristics at time t (of 
the two dimensional system ()2.6|) ) whose value at time was the grid point {xk, yi). The cubic B-spline S 
is defined as follows 

r (2- |:e|)3 if 1 < |x| < 2, 

65(x) = < 4-6x2 + 3|a;|3 ifO<|a;|<l, 
[ otherwise. 

In the expression (|3.8p . the weight is associated to the particle located at the grid point {xk,yi) at time 
i"; it corresponds to the coefficient of the cubic spline determined by the following interpolation conditions 

f{e,x,,y^) = Y.u;liS{x,- X,{r;xk,yi,n)S{y, - X2{r;xt,yi,n) 

k,l 

= ^^kjS{x,- Xk)S{yj -yi). 

k,l 

Adding boundary conditions (for example the value of the normal derivative of / at the boundaries, we 
obtain a set of linear systems in each direction from which the weights ; can be computed as in [191 US] • 
We can now express the full algorithm for the forward semi-Lagrangian method 
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• Step 0: Initialize /^^ = /o(a;j,j/j) 

• Step 1: Compute the cubic splines coefficients wj? ; such that 

/5 = J2'^ti^'^^' ~ ^k)S{yj - yi), 

k,l 

• Step 2: Integrate (|2.6p from to t"+^, given as initial data the grid points = {xk,yi) to get 
X{t;xk,yi,t^) for i G assuming the advection velocity U is known. We shall explain in the 
sequel how it is computed for our typical examples. 

• Step 3: Project on the phase space grid using p.Sp with t = to get f^J'^ = yj) 

• Step 4: Compute the cubic spline coefficients such that 

k,l 

• Go to Step 2 for the next time step. 

3.2 FSL: An explicit solution of the characteristics 

For BSL, especially for the solution of the characteristics, it is possible to choose algorithms based on two 
time steps with field estimations at intermediate times. Generally, you have to use a fixed-point algorithm, a 
Newton- Raphson method (see [IS]), a prediction correction one or also Taylor expansions (see [12]) in order 
to find the foot of the characteristics. This step of the global algorithm costs a lot (see [H]). It is no longer 
needed in FSL, where the starting point of the characteristics is known so that traditional methods to solve 
ODEs, like Runge-Kutta algorithms can be incorporated to achieve high order accuracy in time. Let us show 
the details of this explicit solution of the characteristics, in Vlasov-Poisson and Guiding-Center models. 

In both cases, the dynamical system (|2.6p has to be solved. With FSL, X(t"), U{X{t^),t^) are known. 
You can choose your favorite way of solving this system on each time step, since the initial conditions are 
exphcit. This leads to the knowledge of X(i"+i) and U{X{r+'^),r+'^) so that Step 2 of the previous global 
algorithm is completed. 

As examples of forward solvers for the characteristics curves, the second-order Verlet algorithm, Runge- 
Kutta 2 and Runge-Kutta 4 will be proposed for Vlasov-Poisson, and, as Verlet cannot be applied, only 
Runge-Kutta 2 and 4 will be used for the Guiding-Center model. 

For Vlasov-Poisson, we denote by X{t") = (Xi(i"), X2(i")) = (x", v") the mesh of the phase space, and 
U{X{t^),t") = (w", _B(a;", t")) the advection velocity. The Verlet algorithm can be written 

• Step 1: VA:,/, w^J"^ - w," = ^£;(4!,i"), 
. Step 2: yk,l,x-+'-x- = A^<t^/^ 

• Step 3: compute the electric field at time t"''"^ 

— deposition of the particles a;J!"j"^ on the spatial grid Xi for the density p: p(xi,t"+^) — J2k i^k 
x;j|^), like in a PIC method. 

— solve the Poisson equation on the grid Xi: E{xi,t'^~^^). 
. Step 4: VA:,/,<r-<t^ = f Eixl+\r+'). 

A second or fourth order Runge-Kutta algorithms can also be used to solve the characteristics curves 
of the Vlasov-Poisson system forward in time. The fourth order Runge-Kutta algorithm needs to compute 
intermediate values in time of the density and the electric field. Let us detail the algorithm omitting the 
indices k,l for the sake of simplicity 
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. Step 1: fci = {v^,E{x'\V')) = (fci(l), fci(2)), 

• Step 2: compute the electric field at intermediate time ti: 

— deposition of the particles on the spatial grid Xi for the density p: p{xi,ti) 
{xl + At/2 fci(l))]. 

— solve the Poisson equation on the grid Xii E{xi,ti). 

• Step 3: compute k2 = + ^fci(2), i;(a;" + ^ki{l),ti) 

• Step 4: compute the electric field at intermediate time t2'- 

— deposition of the particles on the spatial grid Xi for the density p: p{xi,t2) = I'^k I'^i^i ~ 
{x]^ + At/2k2{m. 

— solve the Poisson equation on the grid x^: E{xj,t2). 

• Step 5: compute fcs = (u" + ^k2{2),E{x'' + ^fc2(l),t2) 

• Step 6: compute the electric field at intermediate time t^: 

— deposition of the particles on the spatial grid Xi for the density p: p{xi,t3) = J2k I'^k I'^i^i ~ 
(x^ + Atfc3(l))]. 

— solve the Poisson equation on the grid Xji E{xi,t3). 

• Step 7: compute k4 = (w" + At ^3(2), E{x" + At ^3(1), ts) 

• Step 8: - X" = ^ [/ji + 2k2 + 2fc3 + k^] 

In both Verlet and Runge-Kutta algorithms, the value of E at intermediate time steps is needed (step 3 
for Verlet and steps 3, 5 and 7 for Runge-Kutta 4). This is achieved as in PIC algorithms by advancing the 
particles (which coincide at time t" with the mesh in this method) up to the required intermediate time. 
Using a deposition step, the density is computed thanks to cubic splines of coefficients w" on the mesh at the 
right time, and thus the electric field can also be computed at the same time thanks to the Poisson equation. 
Using an interpolation operator, the electric field is then evaluated at the required location (in steps 3, 5 
and 7). Let us remark that this step involves a high order interpolation operator (cubic spline for example) 
which has been proved in our experiments to be more accurate than a linear interpolation (see section 4). 



For the Guiding- Center equation, the explicit Euler method, and also Runge-Kutta type methods (of order 
2, 3 and 4) have been implemented. There is no technical difficulty with computing high order methods. 
This is one of the general interests of forward methods. The time algorithm for solving the characteristics 
at the fourth order is similar to those presented in the Vlasov- Poisson case. However, there is a additional 
difficulty in the deposition step which enables to evaluate the density at intermediate time steps; indeed, the 
deposition is two-dimensional since the unknown does not depend on the velocity variable in this case. 

Let us summarize the main steps of the second order Runge-Kutta method applied to the guiding center 
model of variables X" = (a;", y") and of advection field U{X'\ t") = E-^iX"", T) 

• Step 1: ^"+1 - X" = AtE^iX",^^) 

• Step 2: Compute the electric field at time 

— two-dimensional deposition of the particles on the spatial grid {xj,yi) for the density p: p{xj,yi,t"'~^^) = 
j:k^kS[xj-x-,+']S[y,-y-,+'] 

— solve the two-dimensional Poisson equation on the grid xj: E{xj,yi, tn+i)- 



Step 3: - X" = 



At 



£;-L(x",r) -h£;-L(x"+i,i"+i^ 



Here, the numerical solution of the two-dimensional Poisson's equation is based on Fourier transform 
coupled with finite difference method. See details in Appendix II. 
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3.3 FSL - BSL Cubic Spline Interpolation 

We are going to compare how spline coefficients are computed recurrently, for one dimensional transport 
problems, for the sake of simplicity. 

FSL: deposition principle On our mesh, the grid points Xi = iAx,i = 0,..,Nx at a time n can be 
regarded as particles. Wc have a distribution function which is projected onto a cubic splines basis. Thus, 
we know f{t",x), Vx, then the particles move forward, and we have to compute f{t"'~^^,Xi), i = 0,...,Nx, 
reminding that / is constant along the characteristics, and that the particles follow characteristics between 
r and 

In fact, to each mesh point Xi, a spline coefficient Uk is linked. The thing to understand, is that these 
coefficients are transported up to the deposition phase. The key is then to compute them recurrently as 
follows: 

• Deposition step 

r+\xi) = J2''kSixi-X{e+';xk,n) 

k 

J2 u^S{xi-X{t^+';Xk,n), 

• Update of the splines coefficients w^"*"^ using the interpolation conditions 

1+2 

r+'(^i)= E ^k^'s{xi-xk), 

fe=i-l 

The number of points which actually take part in the new value of /"^^(xi) (here 4) is directly linked with 
the spline degree you choose. A p-Spline for example has a (p + 1) points support. 
In a ID way of regarding the problem, you can easily prove the mass conservation: 

m"+i = AxJ2r^\xi) 

i 

= AxY.Y.'^kSixi- X{r+';xk,n) 

i k 

k i 

Merely with the spline property of unit partition S{x — Xj) = 1 for all x. 

BSL: interpolation principle Let us introduce some notations. The foot of the characteristics Xj, 
belongs to the interval [x; , x^+i [. Then the reconstructed distribution function can be written 

• Interpolation step 

r+^{xi) = r(x(r;xi,i"+i)) 

1+2 

= E ^kS{x{r;xi,r+')-xk) 

k=l-l 

• Update of the splines coefficients to^~^^ using the interpolation conditions 

i+2 

r+\xi)= E ^i-^'s{x,-xk) 

k=i-l 
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Figure 1: Principle of FSL (left) and BSL (right) for linear splines. 

Here, we denoted by X(i"; Xi, f"+^) the foot of the characteristic coming from Xi. The reader is referred to 
[l9l [12] for more details on BSL interpolation. 

In both cases, a linear system has to be solved, of equivalent complexity, so our method is as efficient at 
this level as BSL is. 

3.4 Basic differences between FSL and BSL 

Let us now explain the basic differences between forward and backward semi-Lagrangian methods. In both 
cases, a finite set of mesh points {xm)m=i..N is used, and the values of the function / at the mesh points at 
a given time step are considered. The aim is to find the new values of / on the grid at the next time step 

BSL For BSL, in order to find the (n-|- l)-th value of / at Xm, we follow the characteristic curve which goes 
through Xm, backward in time, until time t". The arrival point will be called the foot of the characteristics 
and does not necessarily coincide with a mesh point. Hence, we use any interpolation technique to compute 
/ at this point, knowing all the values of the mesh at this time. This leads to the new value of f{xm)- Let 
us summarize; 

• find the foot of the characteristics X(t") knowing X(t"+^) — Xm (mesh point) 

• interpolate using the grid function which is known at time t". 

FSL For FSL, the principle is quite different. The characteristics beginning at time t" on the grid points 
are followed, during one time step, and the end of the characteristics {i.e. at time t"+^) is found. At this 
moment, the known value is deposited to the nearest grid points (depending on the chosen method). This 
deposition step is also performed in PIC codes on the spatial grid only, in order to get the sources for the 
computation of the electromagnetic field. Once every grid points has been followed, the new value of / is 
obtained by summing all contributions. The FSL method can be summarized as follows 

• find the end of the characteristics X(t"+^) leaving from X{t'"-) — Xm (mesh point) 

• deposit on the grid and compute the new particle weights. 



4 Numerical results 

This section is devoted to the numerical implementation of the forward semi-Lagrangian method. In par- 
ticular, comparisons with the backward semi-Lagrangian method will be performed to validate the new 
approach. 




RR n° 6727 



10 



Crouseilles, Respaud & Sonnendriicker 



4.1 Hill's equation 

In order to check that high orders are really reached, a particularly easy model can be used, in which there 
are no self-consistent fields. This leads to a ID model with an external force field written —a{t)x, where a 
is a given periodical function. The Vlasov equation becomes: 

df 

-^+vd^f-a{t)xdj = 0, (4.9) 

The solution of this equation is seen through its characteristics, solutions of 

dX dV , , , , 

thus, X is solution of Hill's equation: 

^ + a{t)X^O (4.11) 

Let's note that this equation can be written in a general way ^ = A{t)u, where A is a matrix valued periodic 
function. Since this is a 2D linear system, its solution is a 2D vector space and it is sufficient to find two 
independent solutions. 

Let uj, ip ^ C^(M+,M), with ijj{t) > V< e M+, so that uj is solution of the differential equation 

So u{t) — a;(t)e*'^(*^ and v{t) — uj{t)e~^'^^*''^ are two independent solutions of Hill's equation (see [13] for 
more details) which can be determined numerically. 

For this test case, the initial distribution function will be: 

fo{x,v)^e-i^-^y{x,v) e [-12,12]2. 

The associated solution f{x, v, t) will depend only on ALo{t). In particular, / will have the same periodicity 
as a and lo. This is what will be used for testing the code. For different orders (2 and 4), and different 

Ai, Xrms{t) = J x-^f{x, V, t)dxdv will be displayed on Fig [21 This function should be periodic, and thus 
should reach the same test value Xrms{^) at each period. The error will be measured between the ten first 
computed values and the exact one, for the ten first periods. Then these errors will be summed, so that a 
norm of the error is dealt with: 

fc=10 

err = ^ e^, with = \xrms{'2.kTT) - x„„s(0)|. 

fc=0 

The order of the method is checked in Figure [H Note that = Ny = 1024 to make sure that convergence is 
achieved for the interpolation step. The expected order is achieved for a certain interval. If At becomes 
too small, a kind of saturation happens. This is due to the term in (where h = Ax — Aw) in the 

theoretical estimation of the error for backward methods (|^), which becomes too high and prevents us from 
keeping the correct order. A forthcoming paper will try to do the same kind of error estimation for the 
forward method. 



4.2 Vlasov-Poisson case 

Landau damping The initial condition associated to the scaled Vlasov-Poisson equation ()2.ip - ()2.2p has 

the following form 

/o(a;,'y) = ^Lcxp(-uV2)(l + acos(fca;)), (x, u) e [0, 27r/fc] x R, (4.13) 



INRIA 



Forward semi-Lagrangian Vlasov solvers 



11 




Figure 2: Error as a function of At for RK2 and RK4 (left) and Xrms as a function of time, for At = 27r/25, 
RK2 and RK4 (right). 



where k = 0.5 is the wave number and a = 0.001 is the amphtude of the perturbation, so that hnear regimes 
are considered here. A cartesian mesh is used to represent the phase space with a computational domain 
[0, 27r/A:] x [—VmaxTVmax], Vmax = 6. The number of mesh points in the spatial and velocity directions is 
designated by = 64 and Ny = 64 respectively. Finally, the time step is equal to At = 0.1 and the Verlet 
algorithm is used to compute the characteristics. 

In this context, it is possible to find an exact value of the dominant mode solution of the linearized 
Vlasov-Poisson equation (see Appendix I for some details). The exact electric field corresponding to the 
dominant mode reads 

E{x, t) = 4a X 0.3677e-"-i^33* sin(0.5a;) cos(1.4156t - 0.5326245). 

On Fig. [31 the analytical solution of the norm of the electric field and the implemented one are plotted. 
It can be observed that the two curves are very close to each other. In particular, the damping rate and the 
frequency of the wave are well recovered (7 — —0.1533 and lu — 1.4156) by the method. Similar precision 
is achieved for different values of k leading to different value of the damping rate and of the frequency (see 
Fig. ED. 

The recurrence effect that occurs with the present velocity discretization on a uniform grid, at Tr « 
80 w"^ can also be remarked. This value is in good agreement with the theoretical recurrence time which 
can be predicted in the free-streaming case (see [15]) Tr = 

This test case has also been solved with the "hybrid" method in which the deposition step is only 
performed every T time steps. In all other steps, the remapping (or deposition) step is not performed, 
therefore, it can be linked with a PIC method. As it was already said, it is not really a PIC method since 
the spline coefficients are different on the phase space grid and are updated at each remapping step, whereas 
in classical PIC methods, these coefficients (called weights) are constant equal to J^" where Npart is the 

^ ^ part 

number of particles. On Fig. |4l the electric field is plotted again, for different values of T, and At = 0.1, 
with Nx — Ny — 128 points. As expected, the method works well, even for large values of T. Only a kind of 
saturation can be observed, and it can be seen that values smaller than 2~^^ are not well treated, because 
of the lack of accuracy of the hybrid method. Nevertheless, the results are convincing: the computation gets 
faster as T gets larger, and a good accuracy is still reached. 
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Figure 3: Linear Landau damping for k = 0.5 (left) and k = 0.4 (right) 



Two stream instability This test case simulates two beams with opposit velocities that encounter (see 
Pl I15j ). The corresponding initial condition can be given by 

fo{x,v) — M{v)v'^[l — acos(fca;)], 

with k = 0.5 and a = 0.05. The computational domain is [0, 27r/fc] x [—9, 9] which is sampled by = Ny = 
128 points. The Verlet algorithm is used to compute the characteristics with At = 0.5. 

We are interested in the following diagnostics: the first three modes of the electric field, the electric 
energy l/2||£'(t)|||2 and the time evolution of the phase space distribution function. 

On Fig. m we plot the time history of the first three Fourier modes of the electric field: \Ei\, \E2\, {E^l 
denotes the amplitudes of E{k — 0.5) , E{k = 1) and E{k = 1.5) respectively. We observe that after an 
initial phase, the first mode exponentially increases to reach its maximum at T « 18 0J~^. After this phase 
and until the end of the simulation, a periodic behavior is observed which translates the oscillation of the 
trapped particles in the electric field; in particular, a vortex rotates with a period of about 18 lu~^. The other 
modes \E2\ and |£'3| also grow exponentially and oscillate after the saturation. However, their amplitude 
remains inferior to that of the first mode. Similar observations can be performed for the electric energy 
which reaches its maximum at T « 18u!~^ after an important and fast increase (from t — 8 to t — 18u!~^). 

This test case was also solved with the hybrid method to test the capability in the nonlinear regime. On 
Fig. the first Fourier mode of the electric field is displayed for different T, with At = 0.5, and 128 points 
in each direction. It can be observed that during the first phase, which is a linear one, even for big T, the 
results are quite accurate for all values of T. The hybrid method seems to have more difficulty after this 
linear phase. Numerical noise, one of the main drawbacks of PIC methods can be observed as T gets higher. 
The phenomenon can be understood looking at the distribution function. On Fig. [Hlthe noise clearly appears 
on /. Noisy values quickly reach high values which prevent the method from being accurate enough. They 
are more important in our hybrid method than in classical PIC ones, because of the deposition step, where 
the particles weight play their role. Indeed, if particles with very different weights are located at the same 
place, the deposition does not take into account properly the particles of low weights compared to those of 
heavy weights. On Fig. [SI it can be seen that the vortex which appears at the middle of the distribution and 
should stay there slowly leaves out of the domain. This can be explained as a kind of diffusion. Nevertheless, 
if T remains very little (2 — 4), the results are really good. It can also be noticed that At plays an important 
role, actually, when a smaller At is chosen, the results remain good for bigger T. As an example, if A< ~ 0.1, 
results remain acceptable until T = 16. 
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Figure 5: Two stream instability: Time evolution of the three first modes of the eleetric field (left) and of 
the electric energy (right). 

Bump on tail Next, we can apply the scheme to the bump-on-tail instability test case for which the initial 
condition writes 

fo{x,v) f{v)[l +acos{kx)], 



with 



f{v) = Up exp(— + nt exp 



\v — 



on the interval [0, 207r], with periodic conditions in space. The initial condition /q is a Maxwellian distribution 
function which has a bump on the Maxwell distribution tail; the parameters of this bump are the following 



10(271)1/2'"'' 10(27r)i/ 



= 4.5, ut = 0.5, 



9, At = 0.5. The Runge-Kutta 4 



whereas the numerical parameters are = 128, Ny — 128, Vmax 
algorithm is used to compute the characteristics. 

We are interested in the time evolution of the spatially integrated distribution function 



Fit,v) = 



f{t,x,v)dx, 



and in the time history of the electric energy l/2||i?(t)||^2- For this latter diagnostic, we expect oscillatory 
behavior of period equal to 1.05; moreover, since an instability will be declared, the electric energy has to 
increase up to saturation at i « 20.95 and to converge for large times to 36% of its highest value (see [T51[T7] ). 

On Figures [7] and [51 we plot the electric energy as a function of time. We can observe that oscillations 
appear, the period of which can be evaluated to 1.; then the maximum value is reached at i « 21 and the 
corresponding amplitude is about 9, which is in very good agreement with the results presented in |15j . Then 
the amplitude of the electric energy decreases and presents a slower oscillation due to the particle trapping. 
Finally, it converges to an amplitude of about 2.8 which is very close to the predicted value. However, for 
very large times (at t = 250 w"!), FSL using lower time integrator algorithms (Runge-Kutta 2 or Verlet 
algorithms) leads to bad results (see Fig. [SJ. A very precise computation of the characteristics is required 
for this test and the use of Runge-Kutta 4 is crucial. The results obtained by BSL (see Fig. [7|) are very close 
to those obtained by FSL using Runge-Kutta 4. Moreover, if At increases (up to At w 0.75), BSL presents 
some unstable results whereas FSL remains stable up to At w 1. The use of high order time integrators 
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Figure 6: Two stream instability: Time evolution of the first mode of the electric field (up and left) and 
distribution function at time t = 100 for T = 1,4, 8. 
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Figure 7: Bump on tail instability: time evolution of the electric energy for FSL (left) and BSL (right). 
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Figure 8: Bump on tail instability: time evolution of the electric energy for FSL RK2 (left) and FSL Verlet 
(right). 



leads to an increase of accuracy but also to more stable results as At increases. It has been remarked that 
linear interpolation of the electric field is not sufficient to obtain accurate results with FSL-Runge-Kutta 4 
and cubic spline are used to that purpose. 

Fig. inland Fig. fTUl shows the time development of the spatially integrated distribution function for FSL 
and BSL. We observe that very fast, the bump begins to be merged by the Maxwellian and a plateau is then 
formed at t w 30 - 40Wp ^ 

4.3 Guiding-center case 

Kelvin-Helmoltz instability In order to validate our guiding-center code, we used two test cases intro- 
duced in [T7] and [TT]. 
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velocity velocity 

Figure 9: Bump on tail instability: time development of the spatially integrated distribution function for 
FSL (left) and BSL (right). 




velocity velocity 

Figure 10: Bump on tail instability: time development of the spatially integrated distribution function for 
FSL (left) and BSL (right). 
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Figure 11: Kelvin Helmholtz instability 1: distribution function at time t = 0,30w, 



4.3.1 First test case 

The corresponding initial condition is 

pix,y,t = 0)=p'>iy) + ep\y) 

coupled with Poisson's equation: 

(j) = cf)^{y) + £(j)^{y) coa{kx) 

The instability is created choosing an appropriate which will perturb the solution around the equilibrium 
one (j)'^). Using the the work of Shoucri, we will take: 

y 

p{x, y^t — 0) — sin{y) + 0.015 sin(-) cos(fcx) 

where k — and Lx the length of the domain in the x-direction 
The numerical parameters are: 

Ny^ 128, At = 0.5. 

The domain size has an impact on the solution. The interval [0, 2tt] will be used on the y-direction, and 
respectively Lx = 7 and Lx — 10 ■ This leads to real different configurations: 

With Lx = 7, Shoucri proved that the stable case should be dealt with. That is what was observed with 
this code. 

With Lx = 10, the unstable case is faced. The results prove it on figure Fig. [TT]andfT2l 
For this test case, the evolution of the energy J E^dxdy and enstrophy J p^dxdy will also be plotted on 
Fig. [T31 These should be theoretically invariants of the system. Like for other semi-Lagrangian methods, 
the energy lowers during the first phase, which is the smoothing one, where micro-structures can not be 
solved properly. Nevertheless, the energy is well conserved. Moreover, on Fig. [T31 FSL using second and 
fourth order Runge-Kutta's methods are compared to the BSL method. As observed in the bump-on-tail 
test, the Runge-Kutta 4 leads to more accurate results and is then very close to BSL. However, BSL seems 
to present slightly better behavior compared to FSL-Runge-Kutta 4. But FSL-Runge-Kutta 4 enables to 
simulate such complex problems using higher values of At (see Fig. [H)) . We observed for example that the 
use of Ai = 1 gives rise to very reasonable results since the L^ norm of the electric field E decreases of about 
4%. Let us remark that BSL becomes unstable for At > 0.7. 
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Figure 12: Kelvin Helmholtz instability 1: distribution function at time t = 50, 500 w. 





Figure 13: Kelvin Helmholtz instability 1: time history of norms of E (left) and of p (right). Comparison 
between FSL and BSL. 
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Figure 14: Kelvin Helmholtz instability 1: time history of norms of E (left) and of p (right). Comparison 
of the results for RK4 for different values of At. 

5 Conclusion and perspectives 

In this paper, we introduced the forward semi-Lagrangian method for Vlasov equations. The method has 
been tested on two different models, the one-dimensional Vlasov-Poisson one, and the guiding-center one. 
Different test cases have been simulated, and they are quite satisfying. The results are in some cases a bit 
less accurate, with respect to the conservation of invariants, than with the classical BSL method, but enables 
the use of very large time steps without being unstable and recovering all the expected aims. No iterative 
methods anymore, and high order time schemes can be use in a straightforward manner. The next step will 
be to test the method with the Vlasov-Maxwell model, in which we will try to solve properly the charge 
conservation problem, which is the ultimate goal. We will try to use PIC results about that conservation, 
for example in [I]. We will also try to prove theoretically the convergence of this method. 

6 Appendix I: Linearized Vlasov Poisson and Landau damping 

In classical plasma physics textbooks, only the dispersion relations are computed for the linearized Vlasov- 
Poisson equation. However using the Fourier and Laplace transforms as for the computation of the dispersion 
relation and inverting them, it is straightforward to obtain an exact expression for each mode of the electric 
field (and also the distribution function if needed). Note that each mode corresponds to a zero of the 
dispersion relation. 

The solution of the Landau damping problem is obtained by solving the linearized Vlasov-Poisson 
equation with a perturbation around a Maxwellian equilibrium, which corresponds to the initial condi- 
tion fo{x,v) = (1 -I- ecos(fca;))/\/27re~" Let us introduce the plasma dispersion function Z of Fried and 
Conte ilOj 

_ 2 2 /"' 2 

Z{ri) = \pKe [z — er]i(j])\ where er]i(r]) = — / e* dt. 



We also have, Z'(r]) — —2(r\Z(r\) -I- 1). Then, denoting by E(k, t) the Fourier transform of E and by E{k, uj) 
the Laplace transform of E, the electric field, solution of linearized Vlasov-Poisson satisfies: 



E{k,t) = ^Res^=^^E{k,uj)e 
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where 

, N(k,uj) 
D[k, uj) 

D{k,Lu) = l-^Z'{^), N{k,uj) ^ -^Z(- 



2fc2 'V2fc ' ' ' 2V2fc2 \/2k' 

The dispersion relation corresponds to 

D{k,u) = 0. 

For each fixed k, this equation has different roots ujj, and to which are associated the residues defining 
E{k, t) that can be computed with Maple. These residues that in fact the values 



(6.14) 



Let us denote by lo^ — R&{^j)^ — Ii^i'-^j)^ will be the amplitude of ()6.14p and </? its phase. 

Remark: For each root = + i^i, linked to re"^, there is another: — cj^ + iuji linked to re""''. Then 
keeping only the roots in which uji is the largest, which are the dominating ones after a short time, we get: 

E{k, t) w re''^e~^('^'-+^"*)* + j^e-'''^e-^(-"'-+^"*)* = 2re'^"* cos(w^i - (/?) 

Taking the inverse Fourier transform, we finally get an analytical expression for the dominating mode of the 
electric field, which we use to benchmark our numerical solution: 

E{x, t) « 4ere'^'* sin(fca:) cos{LOrt — 

Remark This is not the exact solution, because we have kept only the highest Laplace mode. Neverthe- 
less, after about one period in time, this is an excellent approximation of E, because the other modes decay 
very fast. 

7 Appendix II: Solution of Poisson in the Guided Center model 

7.1 Find 

We have to solve: 

-A(j>{x,y) = p{x,y) 
We use a Fourier transform in the x direction. This leads, for i G [1, Nx]: 

= ^ <f>iiy) + Piiy) 

Let us introduce a notation. Thanks to Taylor Young formula, we have: 

= T-o = 1 + -TITTT^ -TTT + 



Ay2 y 12 dy^ J dy^ 
Let us apply this to 4>i 

A 2 a2 \ _ _ A 2 

1 + ^TTT + Pi) + 0{Al) = + ft + SrieSc^^ + Spi) + 0{Al). 



12 dyy'^ ^' ' '^'^ ' ' y ^ ^' ' ^' ' 12 
Now, factorizing all of it, we get 



£2A2\ ^ / 10f2A2\ / £2^2 



This is nothing but the solution of a linear system 

where A is a tridiagonal and symmetric matrix and i? is a modified right hand side which allows to achieve 
a fourth order approximation. 
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7.2 Find E 

To compute the electric field from the electric potential, we have to solve E = —Vcp. To achieve this task, a 
quadrature formula is used. 

In the X direction, which is the periodic one, a third order Simpson method is used 

rxi+i 112 

J E{x,y)dx = -(l){xi+i,y) + ^{xi-i,y) w -Ei_i{y) + -Ei+i{y) + -Ei{y) 

where the {4>i)i is given by previous step. There is no problem with extreme values, since the system is 
periodic. We then find the values of the electric field E solving another tridiagonal linear system. 

Whereas on the y-direction, Dirichlet conditions are imposed at the boundary. So, we can use the same 
strategy within the domain, but not on the two boundary points. Since we have a third order solution every- 
where, we want to have the same order there, therefore we cannot be satisfied with a midpoint quadrature 
rule which is of second order. So we will add corrective terms, in order to gain one order accuracy. Here is 
how we do this. 



/ 



VI dv dv^ 

E{x, y)dy = -(f>{x, 1) + ^{x, 0) « -f{E{x, 0) + E{x, 1)) - -^{p{x, 1) - p{x, 0)) + U 
yo ^ 

where 

U = {^e(MD-) - ^))) = % {-dU^ix, 1) - 4>{X, 0))] . 

We want to find the precision of this method, thus, we would like to evaluate the following difference which 
is denoted by A: 

A = cj>{x, 0) - <t>{x, 1) - '^{E{x, 0) + E{x, 1)) + 1) - P{x, 0)) - U- 



Using the Poisson equation and Taylor expansion, we have 



A + E{x,Qi)+E{x,l) 




0(x,O)) 




= 


cP{x,0)) 




2 

dy 


<^(x,0)) 


where e[yo,yi\. Thus, 


, we finally obtain 




A + E{x, 


0)+E{x,l) = -^ 
dy 


{<Pix,l) 



^{dyy{cl>{X,l)-^{X,0)) 

dii'^ 



Moreover, classical quadrature theory gives us the existence of ^2 € [2/0,2/1] such as 

/■^i dv dv^ 92 

Eix, y)dy = -f{E{x, 0) + E{x, 1)) - -^—E{x, 6)- 

Replacing E{x, ^1) by E{x, ^2), which is of first order in our computation leads to 

A + E{x, 0) + E{x, 1) = 1) " <i>{x, 0)) - ^ j"^' E{x, y)dy + {E{x, 0) + E{x, 1)) + ^(A^) 

so that A = O(A^) which is what was expected. 
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